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We exploit the intrinsic difference between disordered and crystalline solids to create systems with 
unusual and exquisitely tuned mechanical properties. To demonstrate the power of this approach, we 
design materials that are either virtually incompressible or completely auxetic. Disordered networks 
can be efficiently driven to these extreme limits by removing a very small fraction of bonds via a 
selected-bond removal procedure that is both simple and experimentally relevant. The procedure 
relies on the nearly complete absence of any correlation between the contributions of an individual 
bond to different elastic moduli. A new principle unique to disordered solids underlies this lack of 
correlation: independence of bond-level response. 


The properties of amorphous solids are essentially and 
qualitatively different from those of simple crystals [1]. In 
a crystal, identical unit cells are interminably and sym¬ 
metrically repeated, ensuring that all cells make iden¬ 
tical contributions to the solid’s global response to an 
external perturbation [2, 3]. Unless a crystal’s unit cell 
is very complicated, all particles or inter-particle bonds 
contribute nearly equally to any global quantity, so that 
each bond plays a similar role in determining the physical 
properties of the solid. For example, removing a bond in 
an ordered array or network decreases the overall elastic 
strength of the system, but in such a way that the resis¬ 
tance to shear and the resistance to compression drop in 
tandem [4] so that their ratio is nearly unaffected. Dis¬ 
ordered materials are not similarly constrained. We will 
show that as a consequence, one can exploit disorder to 
achieve a unique, varied, textured and tunable global re¬ 
sponse. 

A tunable global response is a corollary to a new prin¬ 
ciple that emerges for disordered matter: independence 
of bond-level response. This independence refers not only 
to the dearth of strong correlations between the response 
of different bonds, but also, and more importantly, to 
the response of any specific bond to different external 
perturbations. We will demonstrate this by constructing 
selected-bond-removal networks, where individual bonds, 
or springs, are successively removed to drive the overall 
system into different regimes of behavior, characterized 
by ratios of different mechanical responses. Starting from 
the same initial network, we can remove as few as 2% of 
the bonds to produce a network with a ratio of the shear 
to bulk modulus, G/B, that is either nearly zero (incom¬ 
pressible limit) or nearly infinite (maximally auxetic [5]) 
merely by removing different sets of bonds. Moreover, 
by using different algorithms or starting with different 
configurations, we find that the region within which the 
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bonds are removed can be confined to strips of control¬ 
lable size, ranging from a few bond lengths to the size of 
the entire sample. This has the practical consequence 
that one can achieve precise spatial control in tuning 
properties of the material from region to region within 
the network-as is needed for creating origami [6, 7] or 
kirigami [8] materials. 

We construct networks numerically by starting with a 
configuration of particles produced by a standard jam¬ 
ming algorithm [9, 10]. We place N soft repulsive parti¬ 
cles at random in a box of linear size L and minimize the 
total energy until there is force balance on each parti¬ 
cle. We work in either two or three dimensions and start 
with a packing fraction, </>, that is above the jamming 
density. After minimizing the energy of a configuration, 
we create a network by replacing each pair of interact¬ 
ing particles with an unstretched spring of unit stiffness 
between nodes at the particle centers [11]. We charac¬ 
terize the network by the excess coordination number 
A Z = Z — Z iso , where Z is the average number of bonds 
at each node and Z\ so = 2d — 2 d/N is the minimum for 
a system to maintain rigidity in d dimensions [12]. 

For each network, we use linear response to calculate 
the contribution Bi of each bond i to the bulk modulus, 
B = Bi (see Appendix for details). The distribution 
of Bi in three dimensions is shown in blue in Fig. 1, where 
data are averaged over 500 configurations, each with ap¬ 
proximately 4000 nodes and an initial excess coordina¬ 
tion number AZj n i t j a i ss 0.127 (corresponding to a total 
number of bonds that is about 2% above the minimum 
needed for rigidity). 

Similarly, we can start with the same initial network 
and calculate Gi, the contribution of each bond to the 
angle-averaged shear modulus, G = ]>UG,. (A finite 
system is not completely isotropic, so the shear mod¬ 
ulus varies with direction [13]; we calculate the angle- 
averaged shear modulus, which approaches the isotropic 
shear modulus in the infinite system size limit [14].) The 
resulting distribution for Gi is shown in purple in Fig. 1. 
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Note that the distributions of the bond contributions to 
B and G are continuous, very broad, and non-zero in the 
limit B i: Gi —»• 0. That is, some bonds have nearly zero 
contribution to the bulk or shear modulus while others 
contribute disproportionately. For both B and G, the dis¬ 
tribution decays as a power law at low values of Bi or Gi. 
These power laws are terminated above (£?,) and (Gi) by 
approximately exponential cut-offs. In comparison, the 
distributions for a perfect crystal would be composed of 
discrete delta functions. 



FIG. 1. Bond-level response. Distribution on a log-log scale 
(inset: log-linear scale) of the contribution of each bond to 
the macroscopic bulk and shear moduli, Bi and Gi, for 3 d 
networks with AZj n i t i a i ~ 0.127. Here i indexes bonds. At 
low Bi or Gi, the distributions follow power-laws with ex¬ 
ponents —0.51 and —0.38, respectively. At high values, the 
distributions decay over a range that is broad compared to 
their means, (Bi) and (Gi). 

We next ask if there is a correlation between how an 
individual bond responds to shear and how it responds 
to compression. Do bonds with a large contribution to 
the bulk modulus also have a proportionately large con¬ 
tribution to the shear modulus? Fig. 2a shows the joint 
probability distribution P(Bi, Gi). There is a nearly van¬ 
ishing (but not identically zero) correlation between how 
individual bonds respond to shear and how they respond 
to compression. This is qualitatively different from what 
one would find for a simple crystal. Thus, Fig. 2a il¬ 
lustrates a previously-unrecognized property that is very 
well obeyed by disordered networks: independence of 
bond-level response. 

This new property suggests that one can tailor the be¬ 
havior of the network by selectively removing (pruning) 
those bonds that contribute more or less than the aver¬ 
age to one of the moduli. By so doing, one can decrease 
one modulus with respect to the other. 

First, we consider the known case of rigidity percola¬ 
tion [4, 15, 16], where a bond is picked at random and 
removed. This pruning is repeated until the system be¬ 
comes unstable at A Z = 0. We have implemented a 
slight variation to this procedure: at each step, a bond is 




FIG. 2. Independence of bond-level response. (A) Joint 
probability distribution of Bi and Gi for 3 d networks with 
AZjnitiai « 0.127. There is little apparent correlation be¬ 
tween the response to compression (Bi) and to shear (Gi) 
for a given bond i. (B) The value of G when bonds with 
the largest (purple squares) and smallest (purple circles) Bi 
are removed is nearly indistinguishable from G when bonds 
are removed at random (purple crosses). Similarly, B is very 
similar whether bonds with the largest Gi (blue triangles) are 
removed or bonds are removed at random (blue pluses). 


removed only if each node connected to this bond has at 
least d + 1 remaining bonds in d dimensions. This is the 
condition for local stability of a particle in the original 
jammed packing [17]. As the excess coordination number 
decreases, the bulk and shear moduli vanish together, so 
that G ~ B ~ A Z [4, 15, 16] (see Fig. 2b). Therefore, as 
shown in Fig. 3, the ratio G/B is independent of A Z. 

We now implement the idea of selected -bond removal 
in a variety of ways. First we remove the bond with 
the smallest Bi, namely the weakest contribution to the 
bulk modulus (provided, as above, that each node con¬ 
nected to this bond has at least d + 1 remaining bonds). 
Since the distribution P(Bi) is continuous and nonzero 
as Bi —> 0, the bond removal has almost no effect on 
the bulk modulus. However, since there is little correla¬ 
tion between the contribution of each bond to the bulk 
and shear moduli, there is a much larger effect on the 
shear modulus. The contributions Bi and Gi of the re- 
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FIG. 3. Tuning global response in three dimensions. The ratio 
of shear to bulk modulus, G/B , for four pruning algorithms. 
Error bars (included) are smaller than the symbols. Lines are 
fits to the data over the indicated range and have slopes, from 
top to bottom, of -7.96, -0.01, 1.01, and 1.82. Starting with 
the same initial conditions, we can tune global response by 16 
orders of magnitude by pruning of order 2% of the bonds. 


maining bonds to the moduli are then recalculated and 
the procedure is repeated to remove the bond with the 
smallest Bi. Figure 2b shows that when bonds with the 
smallest Bi are successively removed, the shear modulus 
linearly proportional to A Z. Furthermore, it is quan¬ 
titatively identical, within numerical precision, to when 
bonds are removed at random. The ability to alter the 
scaling of the bulk modulus without affecting the scaling 
of the shear modulus is a clear demonstration that the 
principle of independence of bond-level response allows 
for very precise tuning of global properties. 

Since removing bonds with the smallest Bi has little 
effect on the bulk modulus, we would expect G/B ^ 0 
as A Z —>• 0. As shown in Fig. 3, we find that G/B ~ 
A Z^ B - , with [j,B- = 1.01 ± 0.01. This behavior is iden¬ 
tical to the scaling found in the original jammed sphere 
packings, where A Z is lowered by decompressing the sys¬ 
tem. In decompressing a jammed packing, this suggests 
that the contacts most likely to disappear are those which 
contribute minimally to the bulk modulus, providing the¬ 
oretical insight into why jamming has anomalous G/B 
behavior. 

We can drive the same initial network to the opposite 
limit, G/B —> oo, by successively removing bonds with 



FIG. 4. Tuning global response in two dimensions. The ratio 
of shear to bulk modulus, G/B, for four pruning algorithms. 
Error bars (included) are smaller than the symbols. Lines are 
fits to the data over the indicated range and have slopes, from 
top to bottom, of -5.36, -0.26, 1.27, and 3.05. Starting with 
the same initial conditions, we can tune global response by 17 
orders of magnitude by pruning of order 1% of the bonds. 


the largest contribution to B. As before, independence 
of bond-level response predicts that the shear modulus 
will again decrease linearly with A Z, as we indeed find 
(see Fig. 2b). However, the bulk modulus will decrease 
more quickly, as prescribed by the high B { tail of the 
distribution, suggesting that the ratio G/B should in¬ 
crease. The result of this successive bond-removal algo¬ 
rithm is shown by the blue squares in Fig. 3. We find 
that G/B ~ A Z tJ,B +, where = —7.96 ± 0.01. Thus, 
the increase in G/B occurs with a much steeper power 
law than the decrease of G/B when the bond with the 
smallest contribution to B is removed. This power law 
implies that the distribution P(Bi/ (Bi)) evolves as bond 
pruning proceeds. 

The algorithms mentioned above can be extended in a 
number of ways. For example, one can remove the bond 
with the largest contribution to the shear modulus to 
drive G/B towards zero. In this case, independence of 
bond-level response implies that the bulk modulus would 
respond as if bonds were removed randomly, so that B ~ 
A Z (see Fig. 2b). However, the shear modulus decreases 
more rapidly; we find G/B ~ A Z IJ ' G +, where /tg + = 
1.82 ± 0.01 (purple diamonds in Fig. 3). 

We can also tune two-dimensional networks with equal 
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ease. We construct spring networks in two dimensions 
with approximately 8000 nodes and an initial coordi¬ 
nation number of AZj n j tia i ss 0.047, which is about 
1% above the minimum needed for rigidity. Figure 4 
shows G/B as bonds are pruned towards A Z —> 0 for 
the same four selected-bond removal algorithms as in 
Fig. 3. When bonds with the smallest Bi are removed, 
we find that G/B ~ A Z^ B - with /i#_ = 1-27 ± 0.01. 
This is close to the behavior known for jammed pack¬ 
ings ( G/B ~ A Z 1 ), though it is certainly not as clean 
as in three dimensions. When we prune bonds that 
resist compression the most (largest Bi ), we find that 
G/B ~ A Z flB +, where hb + = —5.36 ± 0.01. At the 
smallest A Z, G/B ~ 10 10 . Finally, when bonds with 
the largest Gi are removed we find that G/B ~ A Z fiG + , 
with /zg + = 3.05±0.01.Although G/B diverges/vanishes 
with slightly different power laws in two dimensions, the 
overall effect is no less drastic. 

Note that our procedures are remarkably efficient in 
tuning G/B. Figures 3 and 4 show that by removing less 
than 2% of the bonds in three-dimensional networks we 
can obtain a difference of more than 16 orders of mag¬ 
nitude in the tuned value of G/B , depending on which 
bonds we prune. In two dimensions, pruning is similarly 
efficient; starting with the same initial configuration we 
are able to obtain differences in G/B that span over 17 or¬ 
ders of magnitude by pruning only ~ 1% of the bonds. We 
also note that our bond-cutting procedures do not cre¬ 
ate any zero-frequency vibrational modes in the system, 
which would herald an instability in the structure. 

The limit G/B —► 0 corresponds to the incompress¬ 
ible limit of a solid where the Poisson ratio, v = (d — 
2 G/B)/[d(d — 1) + 2 G/B] in d dimensions, reaches its 
maximum value of v = +1 (in 2d) or +1/2 (in 3d). The 
limit G/B —► oo corresponds to the auxetic limit where 
the Poisson ratio reaches its minimum value of v = — 1. 
By using these different pruning algorithms, we can tai¬ 
lor networks to have any Poisson ratio between these two 
limits. This ability provides great flexibility in the design 
of network materials. 

We turn now to spatial correlations between cut bonds. 
Driscoll et al. [18] have conducted numerical simula¬ 
tions in which they removed bonds with the largest 
strain under uniaxial or isotropic compression or shear. 
They showed that the cut bonds form a damage zone 
whose width increases and diverges as the initial ex¬ 
cess coordination number, AZi n ;ti a i —» 0; for sufficiently 
small AZj n j t j a i, the pruned bonds are homogeneously dis¬ 
tributed throughout the entire system. Outside this zone, 
they found that the network is essentially unaffected. 

When pruning bonds with the largest contribution to 
B or G, all the data presented thus far are for systems 
with a sufficiently small AZ; n ; t j a i so that the distribution 
of the cut bonds appears homogeneous. In our simu¬ 
lations with large AZi n i t i a i, where the damage zone is 
smaller than the size of our system, we find that G/B 
still diverges/vanishes, but does so when A Z > 0. When 
we remove the bond with the smallest contribution to 


B or G, the bonds are initially removed homogeneously 
throughout the system, independent of AZ; n itiai. The ex¬ 
istence of tunable strong spatial correlations in the cut 
bonds, as found by Driscoll et al. [18], allows one to create 
textured materials spatially varying mechanical proper¬ 
ties. One region may be highly incompressible while a 
nearby region may be highly auxetic. This offers a great 
variety in the mechanical response of these networks. 


For many materials [5] the Poisson ratio decreases with 
increased connectivity of the constituent particles and in¬ 
creases with packing density. We note that neither of 
these correlations hold for the algorithms we have in¬ 
troduced for tuning the Poisson ratio (or ratio of shear 
and bulk moduli). We can reach G/B -+ oo (minimum 
Poisson ratio) or G/B -+ 0 (maximum Poisson ratio) 
by removing the same number of bonds from the same 
starting configuration. Neither the overall connectivity 
nor the overall density is different in the two final states. 
Thus, our procedures for producing tunable Poisson ratio 
materials are fundamentally different from correlations 
considered in the literature. 


We have presented a number of ways of tuning G/B. 
Our results suggest that these ideas may be extended to 
other global properties (e.g., thermal expansion or elec¬ 
trical response [19, 20]) where the response can be writ¬ 
ten in terms of sums over bond contributions. As long as 
there is independence of bond-level response, one should 
be able to tune the ratio of global properties by using 
the same protocol of removing bonds that are especially 
susceptible (or especially unsusceptible) to a given global 
perturbation. 


Our results demonstrate that disordered networks pro¬ 
vide particularly elegant opportunities for constructing 
mechanical metamaterials with tunable, flexible and spa¬ 
tially textured response. However, the algorithms we 
have presented may not be restricted to artificially con¬ 
structed materials. For example, compressing a network 
composed of springs that fail when stressed past a given 
threshold would result in the same network as removing 
springs with the largest Bi, provided that the thresh¬ 
old is sufficiently small. It is also not beyond imagina¬ 
tion that one could selectively break bonds at the nano¬ 
scale level in response to global perturbations in com¬ 
plex solids. Indeed, biology appears to be able to target 
structures in networks that are under particularly high 
stress and to enhance their strength (such as in trabecu¬ 
lar bone [21]). Alternatively, there may be mechanisms to 
buckle or sever strongly stressed fibers (such as in actin 
networks [22]). It is interesting to ask if such selective 
repair or destruction of biological structures changes ra¬ 
tios of different mechanical responses such as the Poisson 
ratio. 
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Appendix A: Calculation of bond-level elastic 
response 

We consider networks of nodes connected by un¬ 
stretched central-force springs with stiffness k = 1. Let 
Sri be the total strain on bond i when the system is de¬ 
formed according to some strain tensor e a p. The change 
in energy of the network is then given to lowest order by 

a s = (ai) 

i 


where Sr .$11 is the component of Sri that is parallel to 
the bond direction. Thus, the bond that contributes the 
most (least) to the response to a given boundary defor¬ 
mation is the one with the largest (smallest) Sr? m. To 
remove the bond that contributes the most to the bulk 
modulus, for example, one would remove the bond with 
the largest Sr? m under compression. This procedure can 
be implemented in either a simulation or an experiment. 

In practice, for our computations, we use linear algebra 
to calculate the response of each bond more efficiently, 
as follows. The bulk elasticity of a system is described to 
linear order by the elastic modulus tensor c a p-ys, so that 
if the system is distorted by the symmetric strain tensor 
e a p, the change in energy is given to leading order by 

A E/V — ^^a/3^a/3 r yS^^5 , (A2) 

where V is the volume of the system. In general, there 
are 6 (21) independent components of the elastic modulus 
tensor in two (three) dimensions, but in the isotropic 
limit this reduces to just the bulk modulus B and the 
shear modulus G. 

The components of c ai a 7 5 are calculated from the 
change in energy of the system under various boundary 
deformations using Eq. Al. The strain Sri can be de¬ 
composed into two distinct parts. First there is an affine 
strain set directly by the strain tensor. However, this 
results in a nonzero net force, f m , on each node to, lead¬ 
ing to a secondary non-affine response. This non-affine 
response is calculated by solving the following system of 
equations 


M. 


i>NA 


= fn 


(A3) 


Under the deformation e a p, the change in energy of bond 
i is 


A Ei — ^ €a0C'i,a(3'y8£ r y5 ■ (A5) 

c i,ap~/8 thus completely describes the bond-level elastic 
response for bond i, and can be used to calculate the 
quantities Bi and Gi considered in the main text. 

The global bulk and shear moduli are linear combina¬ 
tions of the components of the elastic modulus tensor. In 
two dimensions, they are 


B ^ (C-xxxx ~b Cyyyy “b ^C XX yy ) 

G g (^Cxyx-y T Cxxxx “b C yyyy 2 C X xyy ) i 


(A6) 

(A7) 


while in three dimensions they are 


B — g ( C-xxxx T Cyyyy "b C zzzz “b ^Cyyzz "*b 2 C X xzz T 2 C X xyy) 

(A8) 

^ = 15 (3 Cyzyz T 3 C X zxz “b 3 C X yxy 

T C-XXXX ~b Cyyyy T C zzzz Cyy ZZ C-XXZZ C X xyy) • 

(A9) 


Finite disordered systems are never perfectly isotropic, so 
the shear modulus always has some dependence on the 
angle of shear. The above expressions for G represent 
the angle-averaged shear modulus, which reduces to the 
shear modulus in the isotropic limit of infinite system 
size. We calculate the contribution of bond i to the bulk 
and shear moduli in exactly the same way: 


Bi — 4 ( Ci,xxxx “b Ci yyyy T 2 Ci^xxyy) 

Gi — g (4 Ci^xyxy “b Ci jX xxx “b Ci^yyyy 2 Ci’Xxyy') > 


(A10) 

(All) 


in two dimensions, and 

Bi g (.Ci^XXXX “b Ci^yyyy Ci^zZZZ T 2 Ci,yy ZZ T 2 Ci^ X XZZ 4 “ 2 Ci^ XX yy 

(A12) 

Gi — jg (3 Ci^yzyz + 3 Ci^xzxz “b 3 Ci^xyxy 

“b Ci^xxxx “b G^yyyy “b Ci^zzzz C-i^yy Z z C-i,xxzz G,.rxyy) 

(A13) 


in three dimensions. 
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where M. mn is the Hessian matrix and is the non- 
affine displacement of each node. The total strain Si\ 
of bond i is calculated from the sum of the affine and 
non-affine displacements of the two nodes that the bond 
connects. Since A E can be written as a sum over bonds, 
so too can the elastic modulus tensor: 

Cafi'yS — ^ ( G,q/37<5■ (A4) 
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